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Abstract. The "first passage-time" problem is an important problem with a wide 
range of applications in mathematics, physics, biology and finance. Mathematically, 
such a problem can be reduced to estimating the probability of a (stochastic) process 
first to reach a critical level or threshold. While in other areas of applications the 
FPT problem can often be solved analytically, in finance we usually have to resort 
to the application of numerical procedures, in particular when we deal with jump- 
diffusion stochastic processes (JDP). In this paper, we develop a Monte-Carlo-based 
methodology for the solution of the FPT problem in the context of a multivariate 
jump-diffusion stochastic process. The developed methodology is tested by using 
different parameters, the simulation results indicate that the developed methodology 
is much more efficient than the conventional Monte Carlo method, which establishes 
itself as an efficient tool for further practical applications, such as the analysis of 
default correlation and predicting barrier options in finance. 

1. Introduction. In a jump-difFusion process (JDP), the dynamics of underlying 
process have two random components: a continuous diffusion component and a dis- 
continuous jump component [1], in which the jump component can be explained as 
a sudden drop of process's value. The first passage time (FPT) problem for jump- 
difFusion processes has Attracted attention of researchers in such diverse fields as 
queuing networks [2], computer vision [3], target recognition [4]. In the financial 
world, many problems also require the information on the first passage time of 
a stochastic process, for example, in modeling credit risk and valuing defaultable 
securities [1], or in predicting barrier options [8]. Furthermore, it is now gener- 
ally accepted that the geometric Brownian motion model for market behavior may 
produce misleading results in [H [9j , such as mismatching credit spreads on corpo- 
rate bonds, underlying derivative prices. Jump-diffusion processes have established 
themselves as a sound alternative to the geometric Brownian motion model. 

However, if we consider jumps in the process, except for very basic process types 
where closed form solutions are available, when the jump sizes are doubly exponen- 
tial or exponentially distributed [5], or when the jumps can have only nonnegative 
values (assuming that the crossing boundary is below the process starting value) 
[5] , where closed form solutions are available, for most problems we can only resort 
to the numerical procedures. 

Monte Carlo simulation is a very promising candidate in dealing with FPT prob- 
lems. However, in conventional Monte Carlo method, in order to avoid discretization 
bias [7], we need to discretize the time horizon into small enough intervals, and to 
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evaluate the process at each discretized time that is very time-consuming. Recently, 
Atiya and Metwally [HIIS] have developed a fast Monte Carlo- type numerical method 
to solve the FPT problem for jump-diffusion process. 

In many financial problems we have to deal with multiple processes in practice 
and consider their jointly crossing the critical level, so it is very useful to develop 
fast numerical procedure for multivariate jump-diffusion process. In this contri- 
bution, we extend the fast Monte Carlo-type numerical methodology to the more 
general case that covers afhne multivariate processes jump-diffusion. The developed 
methodology can be easily extended to other financial applications and areas where 
FPT problem arises. 

The article is organized as follows: section [2] describes our mathematical model. 
The algorithms are presented in section[3]and simulation results are given in section 
m Concluding remarks are given section [51 

2. Mathematical Model. In this section, first, we give brief discussion on the 
affine jump-diffusion model, then we deduce the first passage time distribution 
under multivariate jump-diffusion process. At last, we consider the issue of kernel 
estimator. 

2.1. Affine jump-diffusion. Affine jump-diffusion is a jump-diffusion process for 
which the drift vector, "instantaneous" covariance matrix and jump intensities all 
have afhne dependence on the state vector. 

Let us consider a complete probability space (SI, F, P) and an information filtra- 
tion {Ft), and suppose that X is a Markov process in some state space D C M", 
solving the stochastic differential equation p]Oj 

dXt = iiiXt)dt + aiXt)dWt + dZt, (1) 

where W is an (i^t )-standard Brownian motion in K"; ji : D ^ M", a : D ^ R"^", 
and Z is a pure jump process whose jumps have a fixed probability distribution v 
on R" and arrive with intensity {\{Xt) : t > 0}, for some A : D — > [0,oo). 

Definition 1. The above model is an afhne model if [10]: 

^l{Xut)^KQ + K^Xt 

{a{Xt,t)a{Xuty\, = {H^\, + {Hi),,X, 

X{Xt) ^ lo + h ■ Xt, (2) 

where K ^ [Kq.Ki) e R" x R"^", H = [Hq.Hi) e R"""" x R«x"X", I = (l^^J^) g 
R"xR"^" and the "jump transform" (determines the jump-size distribution) tpic) = 
/g„ exp(c, z)di'{z), for c G C", is known whenever the integral is well defined. The 
"coefficients" {K,H,l,ip) of X completely determine its distribution. 

2.2. First passage time distribution. Now, we will consider the first passage 
time distribution in the context of multiple processes. In order to obtain a com- 
putable multi-dimensional solutions of FPT distribution, we need to simplify Eq. 
(HI and (Abased on the following assumptions: 

1. Each Wt in Eq. ^ is independent; 

2. Ki ^ 0, Hi — and = that means the drift term, the diffusion process 
(Brownian motion) and the arrival intensity are independent with state vector 
X; 

3. The distribution of jump-size Zt is also independent with X. 
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In this scenario, we can rewrite Eq. ([T]) as 

dXt = iidt + adWt + dZt , (3) 

where 

fj. = Ko, aa^ = Hq, X = ^o- 

Atiya et al. [9] have deduced one-dimensional first passage time distribution in 
time horizon [0,T]. In order to judge multiple processes, from Eq. ([3]), we obtain 
multi-dimensional formulas and simplify them into computable formulas. We will 
highlight the main steps of our procedure below. 

As defined, the multiple processes X can be written a.s X — [Xi, X2, ■ ■ - Y^ ■ Let 
us consider one of its components, a sub-process Xi that satisfies the following 
stochastic differential equation: 

dXi ~ pLidt + aijdWj + dZi 
3 

= fiidt + aidWi + dZi, (4) 
where W, is also a standard Brownian motion and cr, is: 



We assume that in the interval [0,r], times of jumps happen for Xi. Let 
the jump instants be Ti,T2,-- - ,TMi- Let Tq = and Tm^+i ~ T. tj equals 
interjump times, which is Tj — Tj^i. Following the notation of 9J, let Xi{T^) be 
the process value immediately before the jth jump, and Xi{Tj~) be the process 
value immediately after the jth jump. The jump-size is then X,{T+) - X,{T-), 
and we can use this jump-size to generate Xi{T^) sequentially. 

If we define Ai{t) as the event that process crossed the threshold Di{t) for the 
first time in the interval [t, t + dt], we have 

= p{A{t) e dt\x,{T+_,),x,{T-)). (5) 

If we only consider one interval [Tj^i, Tj], we can obtain [TTl [T^] 

Xi(T+)-DAt) 



* exp 



[x,{T;)~-D.,{t)~^H{TJ-t)] 



2 ' 



w^^. ' 



where 



Vi = 7== exp 



[X,(7;+,)-X,(rr)-f;,,r,] 



(Ti yJi'KTj \ 2Tj(T, 



2 



After getting result in one interval, we combine the results to obtain the density 
for the whole interval [0, T] . Let B{s) be a Brownian bridge in the interval , Tj] , 
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with B{T+_^) = X,{T+_^), B{Tp = X.iTp, the probabihty that the mimmum of 
B{si) is always above the boundary level is [13] 



= P M Bis.,)>DMB{T+_,) = X,{T+_,),B{Tr)=X,{T. 



T,-i<s,<Tj 

ifX.(rr) > A(t), 
otherwise. 




2[X.(T+ i)-D,(t)][X,(T-)-D.(t)] 

T":;^ 



Then i?(si) is below the threshold level, which means the default happens or 
already happened, and its probability is 1 — Pij. Then let L{si) = Li denote the 
index of the interjump period in which the time Si falls in [rj,._i, T;,.]. Also let li 
represent the index of the first jump, which happened in simulated jump instant. 

li = min(j : Xi{T^:) > Di{t);k = 1, . . . , j, and 

X,iT+) > AW; fc - 1, . . . , J - 1, and X,{T+) < D,{t)). (8) 

If no such li exists, we set = 0. 

Then we get the probability of the interval [0, T] 

P{A,{si) e ds|X,(T+_i), X,(Tr), j = 1, . . . , + 1) 

giLi {si) Y{k=i Pik if Li < U or li = 0, 

9^LAS^) X{k=l P^k + P^k5{s^ ' W if L, - /„ (9) 

if > 

where 5 is the Dirac's delta function. 

2.3. The kernel estimator. For each X,, after generating series of first passage 
times Si, we use a kernel density estimator with a Gaussian kernel to estimate 
the first passage time density (FPTD) /. As described in [3], the kernel density 
estimator is based on centering a kernel function of a bandwidth as follows: 

^ 1 ^ 

f^^Y.^{h,t-Si), (10) 



N ■ 



where 



The optimal bandwidth in the kernel function K can be calculated as [TS] : 



hopt = (^N^ J in'fdtj , (12) 

where N is the number of generated points and ft is true density. Here we use the 
approximation for the distribution as a gamma distribution as proposed in [5]: 

ft^^f-'^M~at). (13) 



-0.2 



In this case the functional becomes 

5 



/>oo ^ 



W.a,T{2p - i) ^^^^ 



2(2/5-»)(r(/3))2 
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where 

Wi = A^, W2 = 2AB, W3=B^ + 2AC, W4 2BC, W5 = C^, 

and 

A^a^, B = -2a{f3~l), C = (/3 - - 2). 

From Eq. (I14p . apparently, in order to get a nonzero bandwidth, we has con- 
strained (3 to be at least equal to 3. Using this constraint, we can obtain the esti- 

^T^ t 



mates of the parameters a and /3 via the method of moments: t — E{t) — jj X^iLi 



and E{t'^) = ^1 the sample standard deviation is C = y E{t'^) ~ t . The 

estimates are a = and /? = ^ > 3. 

The kernel estimator can be easily generalized to the multivariate case. Sup- 
pose we consider X ~ [Xi, X2, ■ ■ ■ , , let t = [ti, t2, . . . , t^], and st = 
[sii, S2i, . . . , Smi], Sji{j — 1,2, ... ,m) is the first passage time for Xj. Then, the 
multivariate kernel density estimator with kernel K and window width h is defined 

by m 

1 ^ 

where 

K{t) = (2^/i2)-™/2 (^„_LtT? j . (16) 

And if we approximate the true density / as a unit m-variate normal density, 
then the optimal bandwidth hopt is [15] 

l/(m+i) 

(17) 



hopt - ^i/(„+4) 



2m + 1 



3. Algorithms. In section[2l we have built up a multivariate jump-diffusion model 
as describe in Eq. ([3]), and its first passage time distribution was also obtained in 
section [2.21 In this section, we will discuss how to simulate the multivariate jump- 
diffusion process efficiently via Monte Carlo method. In conventional Monte Carlo 
method, the simulation is very straightforward, we divide the time horizon [0, T] 
into n small intervals [0,ti], [ii,t2]5 • ' ■ j \tn-i,T] and in each Monte Carlo run, we 
need to calculate the value of Xi at each discretized time t. We should mention 
that in order to exclude discretization bias, the number n must be large enough. It 
is obvious that this conventional method is very time-consuming. 

Recently, Atiya and Metwally [H |9] have proposed two fast Monte Carlo type 
methods, which are about 10-30 times faster than the conventional Monte Carlo ap- 
proach. We called them uniform sampling (UNIF) method, which involves sampling 
using uniform distribution, and inverse Gaussian density sampling (IG) method, 
which uses inverse Gaussian density method for sampling. 

In this article, we mainly focus on the uniform sampling (UNIF) method and 
extend it to the multivariate jump-diffusion process. The major improvement of 
UNIF method is that it only evaluate Xi at generated jump instants and between 
each two jumps the process is a Brownian bridge, so we just consider the probability 
of Xi crossing the boundary level in {Tj^i,Tj) instead of evaluating Xi at each 
discretized time t. More exactly, we assume that the values of Xi(T^_^) and Xi{T~) 
are known as two end points of Brownian bridge, we generate a variable Si with 
uniform distribution and by using Eq. ([7|) to see whether Xi{si) is smaller than the 
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threshold level, if it defaults, then we have successfully generated a first passage 
time Si and can neglect the other intervals and perform another Monte Carlo run. 

In [51 [3], the jump-diffusion process is involved in a univariate model. However, 
our model is based on multivariate process in which Xi are correlated as described 
in Eq. ([3]), so we need to consider several points as follows: 

1. In this article, as a first step, we assume that the arrival rate A for the Poisson 
jump process and the distribution of (Tj — are the same for each Xi. 
As for jump-size, we should generate them as given distribution, and it can 
be different to reflect the different jump process for each Xi. 

2. At present, we use exponential distribution (mean value fix) for {Tj — T}-i) 
and normal distribution (mean value /xj and standard deviation crj) for the 
jump-size. Of course, we can use any distribution as desired. 

3. If we consider m processes, i.e., X — [Xi,X2, . . . , Xm]^ , then we need an ar- 
ray IsDef ault (whose size is to) to indicate whether process Xi has crossed the 
threshold in this Monte Carlo run. If Xi has crossed, then we set IsDef ault (z) — 
1, and will not evaluate it during this Monte Carlo run. 

Next, we will give a description of our algorithm, based on a multivariate exten- 
sion of the algorithms proposed in [H [9] . 

3.1. Uniform sampling method. Let us consider m processes in the given time 
horizon [0, T], as described above, we have generated the jump instant Tj by gener- 
ating interjump times {Tj — Tj_i), besides we set IsDef ault(z) = 0{i = 1, 2, • • • , m) 
at first. 

From Eq. (U, we can see that, 

1. If jump doesn't occur, the diffusion follows a standard Brownian motion, 
Wi{T) - N{0,T), so interjump size {X,{Tf) ~ X,{T+_^)) follows a normal 
distribution of mean ^.i{Tj — Tj_i) and standard deviation Oi^jTj — Tj_i. 
After extend if necessary, we get 

m 

X,{TT) ^ X,(T+ J + m.(T;- - T,_i) + ^ a.fciV(0, T, - r,_i), 

fe=i 

and the initial state is Xi(0) = Xii^^\ 

2. If jump occurs, the jump-size and direction of Zi{Tj) are not fixed either. We 
simulate the jump-size by a normal distribution, and of course we may gen- 
erate it according to other distribution. Then we can compute the postjump 
value: 

X,{T;)^X,{Tp + Z^,{T,). 

After generating beforejump and postjump value Xi{T~) and Xi{Tp {j = 
1, • • • , M, M is the total number of jumps for all the processes Xi), we can compute 
Pij according to Eq. ([7|). To recur the first passage time density (FPTD) fi{t), we 
need to consider three conditions for each Xi that is still above the threshold: 
1. First passage happens inside the interval. We know if Xi{Tp^) > 
Di{t), Xi{Tj') < Di{t), then the first passage happened in the time interval 
[Tj_i,Tj]. To judge when the first passage happen, first we compute the 
probability Pij of Xi always above the threshold according to Eq. ((7]), then 
we generate Si as Si = bijUi +Tj-i, where bij = -jz^s-^, and u,j is a uniform 
random number in [0, 1]. If Si also belongs to interval Tj], then the first 

passage time occurred in this interval, where s,; is the first passage time and 
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now we set IsDefault(?) = 1 to indicate process Xi has crossed the critical 
level. In this condition, we can get conditional interjump first passage time 
density of that specific interval by Eq. Q . To get the density of whole interval 

[0, r], we have = (^ ^ilp.~^ ) * K{hopt,t - s^), where n is the 

iteration number of Monte Carlo cycle. 

2. First passage doesn't happen in this interval. If Si doesn't belong to 
interval [Tj-i, T^], then the first passage time has not yet occur in this interval. 

3. First passage happens in the right boundary of interval. If Xi{Tj') < 
Di {t) , Xi {T^ ) > Di (t) , which follow the definition in Eq. ^ , then obviously 
Tj. (set li ~ j) is the first passage time. Evaluate the density function using 
kernel function fi,n{t) — K{hopt, t — T].), and set IsDef ault(i) — 1. 

Then we increase j and examine the next interval and judge these three conditions 
for each non-crossing process Xi again. For each Monte Carlo run, if we make a 
rough assumption that the probability of Xi crossing the threshold is not correlated, 

then we can obtain the multivariate FPTD as /„( t ) = 11™ ( "^i-p'^ ) 9ij{^i) * 
K{hopt,~t - ~s). 

After running N times of Monte Carlo cycle, we get the one-dimensional FPTD 
of Xi as fi{t) = ^ Y.n=i fi,n{i)^ and multivariate FPTD a.s f{t) ^ ^ J2n=i /"( * ) 

4. Simulation results. In this section, as a demonstration, we will test the mul- 
tivariate UNIF method on two-dimensional case. In order to check the efficient 
and validity of the UNIF method, we use three examples with different arrival rate 
A = 1, 3, 8 for the Poisson jump process to judge the efficiency of our algorithms. 
The parameters are as follows 

Xq = [0,0]^, D{t) = [ln(0.9) - 0.002t,ln(0.95) - 0.012i]^ 

0.2 0.0 



H = [-0.002,-0.012] 



0.0 0.2 

T 



f^z = [0,0] ' , az = [0.2,0.12] 

where Xq is the starting value for the process, D^^t) is the threshold, /i is the constant 
instantaneous drift, a represents the Brownian motion, and iiz and az are the mean 
and standard deviations, respectively, of the jump sizes. 

The simulation was carried out with total Monte Carlo runs = 500, 000 in 
horizon [0, 1]. Moreover, we have also carried out conventional Monte Carlo simu- 
lation with the same parameters, the estimated density functions are displayed in 
Fig. [T][3l All the simulations were carried out on a 2.4 GHz AMD Opteron(tm) 
Processor. The optimal bandwidth and CPU time are described in Tabled] 

From Fig. [l][3l we can see that the multivariate UNIF method gives similar 
density function as the conventional one, that check the validity of our algorithms. 
An interesting phenomenon is that increasing A will affect FPTD of Xi a lot whose 
probability of crossing the threshold is low in the interval [0, T]. 

Undoubtedly, from Table [l] one can easily realizes that the multivariate UNIF 
approach is much more efficient than the conventional one, which establish it as a 
fast methodology for practical applications. 

5. Conclusion. In summary, we have studied the first passage time problem in the 
context of multivariate jump-diffusion processes. We have extended the fast Monte 
Carlo-type numerical method - the UNIF method - to multiple processes. From our 
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One-dimensional density function 




Figure 1. Example 1 (A = 1): One-dimensional (top) and two- 
dimensional density function (bottom) estimate using 100, 000 it- 
erations for UNIF and conventional Monte Carlo approaches. The 
discretization size of time horizon is A = 0.0002 for conventional 
Monte Carlo method. 



One-dimensional density Function 




Time Time Time " Q Time 

Figure 2. Example 2 (A = 3): One-dimensional (top) and two- 
dimensional density function (bottom) estimate using 100, 000 it- 
erations for UNIF and conventional Monte Carlo approaches. The 
discretization size of time horizon is A = 0.0002 for conventional 
Monte Carlo method. 
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Figure 3. Example 3 (A = 8): Onc-dimcnsional (top) and two- 
dimensional density function (bottom) estimate using 100, 000 it- 
erations for UNIF and conventional Monte Carlo approaches. The 
discretization size of time horizon is A = 0.0002 for conventional 
Monte Carlo method. 

Table 1. The optimal bandwidth hopt, and CPU time per Monte 
Carlo run of the simulations. The first hopt is for Xi and the 
second for X2- All the simulations were performed with Monte 
Carlo runs N — 100, 000, besides, for conventional Monte Carlo 
(CMC) method, the discretization size of time horizon is A = 
0.0002. 







Optimal bandwidth 


CPU time 






Xi 






Example 1 


CMC 


0.012664 


0.006943 


0.286642 




UNIF 


0.016030 


0.013880 


0.000527 


Example 2 


CMC 


0.011157 


0.006582 


0.284554 




UNIF 


0.012249 


0.009443 


0.000731 


Example 3 


CMC 


0.008894 


0.005921 


0.299156 




UNIF 


0.009117 


0.006542 


0.001222 



simulation results, we can see that the multivariate UNIF approach is much more 
efficient than the conventional Monte Carlo method, which illustrates that the de- 
veloped methodology can provide an efficient tool for further practical applications, 
such as the analysis of default correlation and predicting barrier options in finance. 
Acknowledgements. We would like to thank NSERC for its support. 
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